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We propose a concept of a qubit based on two coupled fractional vortices pinned at two artificially 
created « discontinuities of the Josephson phase in a long Josephson junction. Each discontinuity can 
be created by a pair of tiny current injectors with the current 7i n j applied. We map the dynamics of 
the system to the dynamics of a single particle in a double-well potential and calculate the effective 
parameters of this potential. By tuning the discontinuities k oc Ji n j during the experiment we 
are able to control the parameters of the effective double-well potential as well as to initialize the 
fractional vortex molecule to a desired state. The system can be used to study macroscopic quantum 
phenomena involving tailored vortex matter. 

PACS numbers: 74.50.+r, 75.45.+j, 85.25.Cp, 03.65.-w 



I. INTRODUCTION 

In a one-dimensional 0-tt long Josephson junction 
(0-7T LJJ) vortices carrying half of the magnetic flux 
quantum <E>o « 2.07 x 10 _15 Wb appear spontaneously 
at the boundaries between and ir regions. 1 These 
scmifluxons 2,3 have a degenerate ground state of either 
positive or negative polarity. The orientation j" or J, de- 
pends on the direction of supercurrent circulating around 
the 0-7T boundary. The classical behavior of semiflux- 
ons has been studied extensively theoretically 1-5 and 
experimentally. 6-8 

One can also construct "molecules" consisting of two 
or more semifluxons, for example, a molecule with anti- 
ferromagnetically (AFM) arranged scmifluxon and anti- 
scmifluxon situated at a distance a from each other in 
a 0-7T-0 LJJ. Such a molecule exhibits two degenerate 
ground states or 4-t- O nc can switch the molecule be- 
tween f 1 and 44" states by applying a small uniform bias 
current to the LJJ. 9,10 

Such an AFM semifluxon molecule was suggested as 
a semifluxon based qubit . 11 The low energy dynamics of 
the system was reduced to the dynamics of a point-like 
particle in a double- well potential and the parameters of 
this potential and the mass of the particle were related 
to the parameters of the 0-7T-0 LJJ. The parameters for 
which quantum effects emerge were estimated. However, 
in such a qubit the height of the barrier in a double- 
well potential depends only on geometrical parameters 
of the 0-7T-0 LJJ, in particular, on the length a of the 
7r part, which makes it impossible to tune the barrier 
during experiment. 

In a conventional LJJ with the Josephson phase cj> and 
the first Josephson relation j s = j c sin </> one can create 
artificially an arbitrary, electronically tunable k disconti- 
nuity of the phase 4>, so that <f>(x) = n{x) + k H(x) and 
fi(x) is a smooth function without discontinuities and 
H(x) is a Heaviside step function. For n = n this junc- 
tion is equivalent to the 0-7r LJJ described above. Such a 



device was proposed and successfully tested. 12-14 By cre- 
ating discontinuities with k ^ tt (\k\ < 2ir) one can spon- 
taneously form and study vortices carrying an arbitrary 
topological charge of p = — n or p = — n + 27rsgnK that 
automatically appear pinned at the ^-discontinuity. The 
fractional flux associated with a vortex carrying a topo- 
logical charge p is $ = $op/(27r). Fractional fluxons 
can form a variety of ground states, 15 have characteristic 
eigenfrequencies 16 and get depinned by overcritical bias 
currents. 17-19 

In the present paper we investigate an AFM molecule 20 
consisting of two fractional vortices in which the topolog- 
ical charge of the vortices can be tuned electronically by 
changing k. For appropriate parameters of the LJJ we 
obtain an electronically tunable qubit. 

The paper is organized as follows. In Sec. II we de- 
scribe the concept which underlies our idea to observe 
macroscopic quantum effects. In Sec. Ill we introduce the 
model of a LJJ with two discontinuities K\ and n<i sepa- 
rated by a distance a. Based on this model we derive in 
Sec. IV analytical expressions for stationary phases and 
their corresponding energies. In Sec. V we show how the 
state of the system can be manipulated. In particular, we 
discuss how the energy barrier can be reduced to reach 
the quantum regime, which is discussed in Sec. VI. In 
Sec. VII we analyze how the states of the molecule can be 
distinguished. Finally, to keep the paper self-contained 
an appendix provides a collection of relations for elliptic 
functions. 



II. CONCEPT 



Previously 



li 



considered an AFM molecule of two 
scmifluxons in a 0-7T-0 LJJ and found that macroscopic 
quantum effects can be observed if the length a of the ir 
segment is restricted to values 11 a c < a < a c + 0.02. Here 
the scaled length a is measured in units of the Joseph- 
son penetration depth Aj. The critical distance a c is a 
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bifurcation point: for a < a c there exists only one static 
solution with constant phase (flat phase state), while for 
a > a c there are two degenerate solutions (AFM states) 
that can be used as two classical states of a qubit. In 
the present paper we only consider infinitely long LJJs, 
where a c = n/2. 

It is rather difficult to fabricate the junction with the 
value of a that lays within the narrow range mentioned 
above. Therefore, it is highly desirable to make a or 
a c tunable. The scaled length a is fixed by design and 
can vary only slightly due to variations of Xj, e.g., with 
temperature. The value of a c , however, can be varied in 
a larger range as described in the following. 

Consider the 0-re-2re LJJ of Fig. 1(a) with two discon- 
tinuities K, see Fig. 1(b). For k = tt this reduces to the 
previous case of a 0-ir-2ir LJJ or, as the Josephson phase 
is 2ir periodic, to a 0-7T-0 LJJ. For a > a c the two classi- 
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FIG. 1. (Color online) (a) Sketch of a 0-re-2re LJJ with a = 2, 
where the two discontinuities (b) are located at x — ±a/2. 
The phases fi(x) (c) and corresponding magnetic field profiles 
Hx{x) (d) describe two states (re, re — 2-7r) and (re — 2n, re). The 
coordinate x is normalized to Aj. 



cally degenerate ground states arc molecules consisting of 
vortices with topological charges (+tt, — it) or (— w, +ir). 
Now, let us vary k around k = tt. The two corresponding 
states shown in Figs. 1(c) and (d) then are (k, k — 2ir) 
and (k — 2tt, k) 15 . Note that the critical distance 15 



(1) 



where F(<f>,k) is the elliptic integral defined in (A.l), 
which is not constant but depends (weakly) on k as shown 
in Fig. 2. We find that the value of a c may change from 
a c (n) = tt/2 « 1.57 to a c (0) = a c (27r) = 21n(l + w 
1.76. 

Thus, we can fabricate a LJJ with the scaled distance 
1.57 < a < 1.76 between the discontinuities. Then, we 
can decrease n starting from k = tt to bring o c (k) as close 
as needed towards a, as shown in Fig. 2. Note that the 
required precision for a is less restrictive by one order of 
magnitude than in the previous proposal 11 . In addition 
we obtain a system which can be tuned into the quantum 
regime. In the following sections we analyze this system 
in detail. 




FIG. 2. The critical distance a c for asymmetric AFM 
molecules as a function of re. It separates the flat phase state 
(white region) from the AFM state (gray region). As an ex- 
ample we consider a LJJ with a = 1.7 (straight line). By 
decreasing re we can tune a c towards a, e.g., re = 0.37T cor- 
responding to a c = 1.66 as indicated by the point on the 
straight line. 



III. MODEL 

The dynamics of the Josephson phase n(x,t) in in- 
finitely long LJJs without bias current and dissipation is 
described by the sine-Gordon equation 3 

Hxx{x, t) - Hu(x, t) - sin[/i(x, t) + 0(x)} = 0. (2) 

We use dimensionless quantities, where lengths are writ- 
ten in units of the Josephson penetration depth Aj and 
time is written in units of u>~ 1 , where uj p is the plasma fre- 
quency. Energies are measured in units of EjXj, where 
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Ej is the Joscphson energy per length. The subscripts x 
and t denote partial derivatives with respect to coordi- 
nate and time, accordingly. 
The function 

e(ar) = -«iH(ar+|)-KaH(aj-|) (3) 

describes two discontinuities —K\ and —k.2 located at x = 
—a/2 and x = +a/2, where H(x) is the Heaviside step 
function. The values of «i and Ki can be controlled 13 
by injector currents Ii n ji oc K\ and Ii nj 2 oc K2- The two 
discontinuities divide the x-axis into three parts to which 
we will refer to as the left, middle and right part of the 
junction, see Fig. 1(a). 

The boundary conditions for /i(x, t) are given by 

H x (-oo,t) = fM x (+oo,t) = 0. (4) 

Although there are discontinuities at x = ±a/2 the phase 
/J,{x,t) and its derivative /j, x (x,t) with respect to x have 
to be continuous at x = ±a/2 at any time t. 

The sine-Gordon equation (2) can be derived from the 
Lagrangian density 

C=y 2 t {x,t)-U, (5) 

where 

U = \nl{x, t) + 1 - cos \p(x, t) + B(x)] (6) 

is the potential energy density. 

To avoid infinite energies the restriction 

cos[/A(±oo,t) +0(±oo,t)] = 1 (7) 

assures that the potential energy density (6) vanishes at 
x = ±oo. This condition restricts the phases at x = ±oo 
to the values 

A*£ = n(—oo,t) = 2ttu l (8a) 
/iR = ^i(+oo, t) = 2nn R + Ki + n 2 , (8b) 

where and hr are integer numbers. 

IV. STATIONARY SOLUTIONS 

In this section we first derive analytical expressions for 
the stationary Joscphson phase (J.(x) and the correspond- 
ing energies. These expressions depend on the values 

Hi = n (-a/2) (9a) 

and 

/i. r = (i i+a/2) (9b) 

of /i at the boundaries (the left and the right) between 
the regions which have to be determined sclfconsistcntly. 
Furthermore, we analyze the stability of the stationary 
solutions. 



A. Phases 

In order to find stationary solutions of Eq. (2) we in- 
tegrate its static version 

fi xx - sm(fi + 6)=Q (10) 

in each region of the junction. This procedure leads to a 
conservation relation 

i/4 + cosG" + 0) =2fc 2 -l, (11) 

or 

where k is a positive real number and has different values 
in the inner and outer parts of the junction. 

By comparing this equation with Eq. (A. 6) we see that 
the phase 

fi(x) = (2n - 1)tt - 6(x) ± 2 am [k(x + c), fc" 1 ] (13) 

solves Eq. (10) in every region. Here n is an integer 
number and c is an arbitrary constant. Furthermore we 
have introduced the Jacobi amplitude function am(u, k) 
defined by Eqs. (A. 4) and (A. 5). 

From the boundary conditions (4) and (7) in Eq. (11) 
we find 

k = 1 (14) 

in the left and right part of the junction. 

By matching Eq. (11) for the outer part (k = 1) and 
the inner part {k ^ 1) of the junction at x = ±a/2 we 
obtain 

,2 i . • K i • / K i\ 
k = 1 + sm — sin - — J 

= l-siny-sin (/"r - Ki - -y)' ( 15 ) 

Therefore, k can vary between the values and \[2 in 
the middle of the junction. Since the behavior of the 
Jacobi amplitude function depends on the value of k as 
discussed in the appendix, we consider three different 
cases: (a) outer regions with k = 1, (b) inner region with 
k > 1 and (c) inner region with k < 1. 

a. Outer regions with k = 1. With the help of 
Eq. (A. 7) we obtain from Eq. (13) the expression 

»{x) - Hl,r + 4 arctan (c a/2 -^ tan W ' r ~ ML ' a )(16) 
for k — > 1, where \ii and /z r are restricted to 

lw,r - Hl.r\ < 27T. (17) 

Since we have to match these two solutions and their 
derivatives to their counterparts in the inner region we 
define the signs 

cr; = sign[/i x (-a/2)] = sign(^/ - // L ), (18a) 
a r = sign[/i x (+a/2)] = sign(^ fl - /i r ) (18b) 

of n x (x). 
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b. Inner region with k > 1. With the help of 
Eq. (A. 4) it can easily be verified that 

fi(x) = Ki — it + 

2a am [k (x + a/2) + aF(ipi, fc" 1 ), fc" 1 ] (19) 

matches the solution for the left part of the junction, 
that is fi(—a/2) = \ii and a = a\. Here we have used the 
abbreviation 



i>l,r = h (Pl,r +7T- Kl) 



(20) 



To match the solution for the right part of the junction, 
that is //(a/2) = fi r and a = a r , we have the additional 
constraint 



a = ak- x [F(Vy,fe _1 )- F(V>z,fc -1 )] . 



(21) 



c. Inner region with k < 1 . Similar to the case k > 1 
we find 

//(#) = (277, — 1)7T + Kj + 

2cr ; am[fc(a; + a/2) + a t fc F(<^, fc), (22) 
where we have used 

ipi^r = arcsin (fc" 1 sin-0;^) (23) 

and 

^j, r + 777T = ^, r , |^ >r | < tt/2. (24) 

Note that there is only one integer number 77. Therefore, 
the values of ipi,r and fj,^ r are restricted to \ip r < ^ 
and |// r — fj,i\ — ^ 7r - 

The additional constraint to match the solution for the 
right part reads 



[a r F(ip r ,k) -aiF(tpi,k)] 
+ (1 - aio- r )K(k) +AvK{k). 



(25) 



Due to the properties of elliptic integrals, for a given 
value of a the integer number v is limited to the values 



a 1 
°^ < 2^ + 2- 



(26) 



We now have expressed the stationary solutions of the 
sine-Gordon equation (2) in terms of //£, fin, Hi and 
fi r . The possible values of //£ and //fl are restricted by 
Eqs. (8) and depend on the topological charge of the sys- 
tem. The values of /// and fi r have to be determined from 
Eqs. (15), (21) and (25). Examples will be presented in 
Sec. VB. 



B. Energies 

In this section we use the results of the previous section 
to calculate the energy of the stationary solutions. 
Integrating Eq. (6) and using (11) we obtain the energy 



Ui = 



-a/2 



fx\ dx 



fi x dfi 



(27) 



of the solution in the left part of the junction. From 
Eq. (16) we find 

/z x = 2sin[i(/z(a;)-/Zi)] (28) 

and arrive at 

Ui = 8 sin 2 [|( w -fi L )} . (29) 

Similarly, we obtain the expression 

U r = 8 sin 2 [\{fi r -fi R )] (30) 

for the right part of the junction. 

For the middle part of the junction we again combine 
Eqs. (6) and (11) and obtain the energy 



[ a 2 [fi 2 x + 2(1 - k 2 )] dx 

J -a/2 



(31) 



of the phase in this region of the junction. From Eq. (13) 
we find 



& = 4fc 2 dn 2 [kix + c),^ 1 ] 



(32) 



and 



± sn [k(x + c), k- 1 ] = (-1)™ cos [| (n(x) - . (33) 

Equations (32) and (A. 11) allow us to evaluate the 
integral. Finally, we obtain with the help of Eqs. (A. 15) 
and (33) the energy 



U m = 2o(l - k 2 ) + 4k£ (ka, fc" 1 ) 

2 2 



4fc 1 cos ^— : — - cos f ~' — - sn(fca, k 1 ),(34) 



where £ is the Jacobi cpsilon function defined by 
Eq. (A. 11). The total energy 



U=U t + U m + U r 



(35) 



of a stationary solution fi(x) is the sum of the energies 
of these three domains. 



C. Stability analysis 

To study the stability of a stationary solution fi(x) we 
insert the ansatz 



/t(x, t) — fi(x) + ip(x)e 



(36) 



into the sine-Gordon equation (2) and linearize it around 
the stationary solution fi[x). Since fx(x) solves the sta- 
tionary sine-Gordon equation (10), we obtain the differ- 
ential equation 

- ip"(x) + cos \p,(x) + 6{x)] i}>{x) = uj 2 ij{x) (37) 

for the eigenmodes ijj(x) with eigenvalues w 2 , which, after 
substituting Eq. (13), takes the form 

i/)"(x) + [l+u 2 -2 sn 2 (fc(x + c), fc -1 )]^^) = 0. (38) 
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Note that Eq. (38) is the Lame equation 21 . The boundary 
conditions and the matching conditions at x — ±a/2 for 
ip(x) follow from the boundary and matching conditions 
for fj,(x) discussed in Sec. III. 

As long as the lowest eigenvalue Wq of the Lame equa- 
tion (38) is positive, the stationary solution ti(x) is sta- 
ble. When it becomes negative, the stationary solution 
H(x) is unstable. In Sec. VB we solve Eq. (38) numeri- 
cally. 

V. TAILORING MOLECULE STATES 

Our ultimate goal is to have a tunable qubit based 
on an AFM molecule. For this purpose we first bring 
the system into one of its degenerate ground states and 
choose k such that the barrier between the two states is 
high. By tuning k we then reduce this energy barrier to 
reach the quantum regime. 



Starting from the k = state and no vortices in the J J, 
one slowly increases the value of k and finds the sym- 
metric ferromagnetic vortex configurations (k, k). When 
K is increased above n])(a), 15 where tt < k]) (a) < 2tt, 
the molecule emits one fluxon and reconfigures into one 
of the degenerate asymmetric antifcrromagnctic states 
(k, k — 2tt) or (k — 2tt,k). To reach the (+7T, — tt) or 
(— 7r, +7r) state we decrease k to tt. Then we further de- 
crease k to reach the quantum limit. This technique has 
two advantages: (i) one needs only one injector current 
source; (ii) one can operate the qubit at small k values, 
which brings less noise. This technique, however, has 
two disadvantages: (i) it does not work in annular LJJs, 
where the emitted fluxon cannot leave easily; 22 (ii) we 
do not know if we prepared the (k, k — 2tt) or (k — 2tt, k) 
state. 

For our further calculations we choose the second ap- 
proach where n, and therefore, the additional noise in the 
system is small in the quantum limit. 



A. Initializing the AFM state 

In a previous publication 11 it was shown that a con- 
figuration with the discontinuities «i = tt and K2 = 
(0-7T-0 LJJ) provides an effective double-well potential 
with degenerate ground states where quantum tunneling 
can be observed. 

An intuitive way to arrive at this state from the 
k\ = K2 = state is to use 

ki = Hi, n 2 = —Ki (39) 

and sweeping Ki from to tt. Then, we want to deal with 
a molecule consisting of a direct and a complementary 
vortex. Therefore, we have to change the discontinuities 
appropriately, e.g., 

Kl = TT + K t , K2 = -TT + K t , (40) 

where the (de)tuning Kt starts at and changes in the 
range between —tt and +tt. 

In this approach one can employ two injector current 
sources: one creating Ki wired to produce two disconti- 
nuities of opposite signs like in Eq. (39) and another one 
creating Kt and wired to produce two discontinuities of 
the same sign like in Eq. (40) . The advantages of this ap- 
proach are: (i) instead of using the first current source, 
one could replace the junction by a conventional 0-7T-0 
LJJ and then only apply injector currents for the tuning 
Kt according to Eq. (40); (ii) this technique works also in 
an annular JJ; (iii) after Ki has reached the value tt we 
know that we have prepared the (+tt, —tt) state rather 
than the (—tt, +tt) state. The disadvantage of this ap- 
proach is that one always injects a current oc k = tt or 
more into the system, which brings additional noise. 

An alternative way is to use only one injector current 
source, which is coupled to injectors as 

Kl = K2 — k. (41) 



B. Tuning the barrier height 

We consider the particular injector-tuning technique 
described by Eq. (41) for characterizing the states, its 
energies and stabilities during this process. 

1. Determination of phase parameters 

After the initialization procedure discussed in the 
previous section the two discontinuities have the value 
k = tt and the AFM state has the topological charge 
fJ-R ~ Ml — 0. For simplicity we choose nj, = and 
rift = —1 in Eqs. (8). 

When we tune k no (anti)fluxon is emitted or absorbed. 
Therefore, til and ur do not change and the topological 
charge 

fi R - hl = 2k-2tt (42) 

is induced in the system when the discontinuities have 
the value k. 

In the stationary solutions presented in Sec. IV we 
still have to find and ii r . In our examples we restrict 
ourselves to the values of a covering the gray region in 
Fig. 2 near the boundary a c («), that is, 1.5 < a < 2 and 
< k < 2tt. It turns out that for these parameters there 
are only solutions for k < 1. For this limitation, Eq. (15) 
restricts the values of Hi and \i r to the intervals 

^ - TT </i;+27T77, 1 < ^ (43a) 

and 

y <fl r +2TTn 2 < tt + y . (43b) 

Here we use instead of Eq. (17) the more restrictive con- 
dition \ni, r — ^l,r\ < Ti" which requires ri\ — and n-2 = l. 
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This condition is reasonable because even in the case 
when we have a whole integer fluxon at one of the discon- 
tinuities the phase changes at most by it when one goes 
from infinity up to the fluxon center. 

By using the constraints Eqs. (43) together with 
Eq. (15) we express /x r in terms of \xi and obtain 



Mr 



-m 



2re- 2?r 



and 



Hr = fll + K — 7T 



(44a) 



(44b) 



as the two only possible solutions. 

Next we substitute these two expressions into Eq. (25) 
and solve the resulting equations numerically for /i; for a 
given distance a. We denote the solution corresponding 
to Eq. (44a) by an index m. and the two solutions cor- 
responding to Eq. (44b) by an index + and — . Figure 3 
illustrates this process for a = 2. 

When we insert the three values of \x m and /z± into 
Eqs. (16) and (22) we find the phases jJ, m {x) and jJ,±(x) 




FIG. 3. (Color online) Possible values of [i\ for a = 2 and 
k = 0.77T. When we insert Eqs. (44a) and (44b) into Eq. (25) 
we obtain for v = two functions ao(fii) depicted by a dotted 
line and by a dashed line, respectively. For a given value a = 2 
three values for fii are possible: m = fj, m for Eq. (44a) and 
fit = /u± for Eq. (44b). 
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FIG. 4. (Color online) The phases h L m {x), fi+(x) and fJ,-(x) 
corresponding to the values fii = fi m , fii — (i+ and fii = fi- 
(arrows in Fig. 3) for a = 2 and k = 0.7-7T. 



depicted in Fig. 4. These solutions satisfy the symmetry 
relations 



and 



H m (—x) = —fi m (x) +2k — 2% (45a) 



H±(—x) = —fJ,^(x) + 2k - 2ir. (45b) 



The solution fi + (x) describes a molecule consisting of a 
k vortex at x = —a/2 and the complementary k — 2tt vor- 
tex at x = +a/2. The solution fJ.— (x) is a complementary 
molecule (re — 2-7T, re), which has the same energy. These 
two solutions are stable whereas the solution fi m {x) is 
unstable, see Sec. VB2. 



2. Energies and stability analysis 

Depending on the parameters of the L JJ there is either 
a single stable stationary solution /i m (x), or two stable 
solutions /J,±{x) and one unstable solution /j, m (x). The 
two stable solutions are separated by an energy barrier 
which is governed by the unstable solution. To reach the 
quantum limit we want to reduce this energy barrier by 
tuning k. Therefore, we investigate how the energies and 
stabilities of the stationary solutions depend on re. 

For our examples we use the following JJ parameters: 



UJ,, 



E, 



2ir 



27r/i 2A L j c 



(46) 

where Ax, is the London penetration depth, j c is the crit- 
ical current density of the J J, C is the capacitance of 
the J J per unit of area, w is the JJ's width, and $o ~ 
2.07 x lCT 15 Wb. Typical parameters are Xl = 90 nm, 
w = 1/irn, and C = 4.1^tF/cm 2 with j c = 100 A/cm 2 . 
This corresponds to \,j — 38 fixa, oj p — 2ir x 42.8 GHz, 
and EjXj = 78.4 mcV. 

Without losing generality we restrict ourselves to 
< re < 7r. The results for the range n < k < 2ir are 
obtained by substituting re — > 2tt — re. 

In order to calculate the energies of the solutions we 
insert the values for [ii obtained in the previous subsec- 
tion into Eq. (35) and obtain the energies U m (n) of the 
state n m (x) and U + (k) — L/_(re) of the states H±(x). In 
Fig. 5 we depict the corresponding energy difference 



AU(k) = U m {n) - U±(k) 



(47) 



for different values of a. This figure shows that the en- 
ergy barrier can be tuned down to zero by decreasing the 
discontinuity re from n to 0. 

The stability of the solutions /J, m (x) and (jl±(x) is deter- 
mined by the eigenvalues of Eq. (38). These eigenvalues 
are calculated numerically. We have used a junction of 
length I = 20 to emulate an infinitely long JJ. The results 
arc depicted in Fig. 6. Positive eigenvalues characterize 
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stable solutions, whereas negative eigenvalues character- 
ize unstable solutions. 

In the lower part of Fig. 2, that is for a < a c (ir) w 1.57, 
there is only one stationary solution fj, m (x), so that 
AU(k) = in Fig. 5. Figure 6 shows that the corre- 




k/tt 



FIG. 5. (Color online) Energy difference AU{k), Eq. (47), 
for different values of a. By tuning k, AU(k) can be made as 
small as needed. The right vertical axis shows the correspond- 
ing temperature AT = AU-EjXj /ftB for typical experimental 
parameters mentioned after Eq. (46). 
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FIG. 6. (Color online) Lowest eigenvalue ujq(k) for the sta- 
tionary solutions fi m {x) and n±(x). When luq(k) is positive, 
the stationary solution is stable; when ujq(k) is negative, it is 
unstable. 



sponding eigenvalue Wq is positive, therefore this solution 
is stable. 

For a c (ir) < a < a c (0) ~ 1.76 we have to distinguish 
between k > K c (a) and k < ft c (a), where the function 
K c (a) is the inverse of a c (ft) defined by Eq. (1). For 
k > n c (a) there are three stationary solutions jJL m {x) and 
H±(x). As shown in Fig. 5 the energy barrier AU(k) 
becomes lower when we decrease k and vanishes at k = 
ft c (a). Our numerics shows that AU(k) rx (k — k c ) 2 for 
k — y k c . From the eigenvalues Wo(ft) displayed in Fig. 6 
corresponding to these stationary solutions we conclude 
that n m (x) is unstable whereas ^±{x) are stable. 

At k = K c (a) all three solutions have the same energy 
and AU(k) vanishes. The eigenvalues 0Jq(k) join at the 
bifurcation point k = n c (a) and vanish, see Fig. 6. For 
k < K c (a) the two stable solutions n±(x) disappear, while 
the unstable solution /x m (x) becomes stable, see Fig. 6. 

Finally, for a > a c (0) there are always three stationary 
solutions n, n (x) and n±(x). According to Fig. 6, n m {x) is 
unstable and n±(x) are stable. All three solutions reach 
the same energy at k — > and the energy barrier van- 
ishes as AU(k) oc «, see Fig. 5. The eigenvalues ujq(k) 
corresponding to the three solutions vanish too at k — s- 0. 

In fact, all three solutions are just single fluxons. The 
two stable solutions fJ,+ (x) and fJ—(x) are fluxons weakly 
pinned at —a/2 and +a/2, respectively. The unstable 
solution /i m (x) is a fluxon at x = 0. It is rather interest- 
ing that AU(k) vanishes for k — > even for a > a c (0). 
Therefore, one can always make the barrier as small as 
needed for arbitrary a > a c (7r) = n/2 and is not limited 
by the interval of 1.57 < a < 1.76! 

This limit k — > is similar to the heart-shaped qubit, 
in which a fluxon is trapped in a double-well potential 
created by a non- uniform magnetic field 23 . In this situ- 
ation, residual pinning in the system can play a major 
role and make the parasitic pinning potential larger than 
the one we construct here. 



VI. QUANTUM REGIME 

We expect to observe coherent quantum oscillations 
between the two degenerate states corresponding to the 
classical states /z+(x) and H-{x) in the region where the 
energy barrier in Fig. 5 becomes sufficiently small, i.e., 
the coupling between these two states is sufficiently large. 
In order to quantify this effect we map the dynamics 
of the system to the dynamics of a single particle in a 
double-well potential and calculate the energy splitting 
5e separating the two lowest energy levels. Because the 
oscillation frequency between the two states n+(x) and 
fi-(x) is given by 5e/h we are then able to determine the 
values of the discontinuity k oc /j n j for a given value of a 
where the quantum oscillations become observable. 
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FIG. 7. The energy differences feoi = £1 — £o (solid lines) and &12 = £2 — £1 (dashed lines) calculated from Eq. (53) and the 
semiclassical expression A (gray lines), Eq. (55), as a function of re for different values of a. Similar as in Fig. 5 the right vertical 
axis shows the corresponding temperature for typical parameters mentioned after Eq. (46). For the parameters indicated by 
two arrows more details are shown in Fig. 8. 



A. Single-mode approximation 

In order to map the dynamics of the complete system 
to the dynamics of a single particle in a double-well po- 
tential we express the phase 



H(x,t) = n m {x) + ^2 q n (t)ip n (x) 



(48) 



in terms of the stationary solution // m (x) and the cor- 
responding eigenmodes ipnix) defined by Eq. (38). By 
inserting this expansion into the Lagrangian density (5) 
and integrating over x we obtain a Lagrangian for the 
mode amplitudes q n (t) which describes the motion of a 
fictitious particle in many dimensions. 

If the lowest eigenvalue w 2 of Eq. (38) is sufficiently 
separated from the next higher eigenvalue we can expect 
that for low energies the particle only moves along the 
"<7o direction" . Motivated by this simple picture, we omit 
the higher modes in Eq. (48) and use the approximation 



fi(x, t) ss n m {x) + q {t)ip (x), 



(49) 



where ipo{x) is the eigcnfunction corresponding to the 
lowest eigenvalue Wg. We insert this equation into our 
Lagrangian density (5) and take into account only terms 
up to the fourth order in q and obtain the Lagrangian 



where 



L = 



U(q ) 



y 2 - U(q ), 



2 2 



21 



(50) 



(51) 



is the potential energy Here the positive parameter K is 
defined by 



K 



+ 00 



dx cos[/i (z) + B{x)] ipofa) (52) 



and ipo is normalized according to f tpQ dx = 1 . Due to 
the symmetry relation (45a) for /j, m (x) there is no third- 
order term in Eq. (51). 

Note that for ljq < 0, where /x m (x) is unstable, U(qo) 
describes a double-well potential with a maximum at 
<7o = and two minima at qo = K\loq\. For 

cjq > it only has one minimum at go = 0. In the 
first case the oscillation frequencies around the minima 
are w± = a/2 \ujq\- In the second case ujq is the oscillation 
frequency around the mimimum. 



B. Energy splitting 

For the Lagrangian (50) we can derive the stationary 
Schrodinger equation 



"off ° 

2 dql 



U(q ) u-j(qo) = £jUj(q ), 



where 



tujjp 



(53) 



(54) 



is the dimensionlcss Planck constant and the energy 
eigenvalues Ej are given in units of EjXj. 

In order to calculate the energy splitting we solve 
Eq. (53) numerically for the potential U(qo) given by 
Eq. (51). Additionally, we compare our numerical results 
for Seqi = £i — Sq with the scaled energy splitting 



8S e ff4 



/ 2AU 
irh c gui± 



exp 



WAU 

3ft G ffW± 



(55) 



in the semiclassical limit. 24 

Figure 7 shows our numerical results for fe i and the 
semi-classical expression A for different values of a. We 
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note that the semi-classical expression A is a good ap- 
proximation for Seqi as long as Seqi is small. To have 
a good quantum-mechanical two-level system (qubit) the 
two lowest eigenvalues £o and £1 have to be well separated 
from the higher eigenvalues. Therefore, we additionally 
compare Sei2 = e 2 — £1 to Seoi- 

To establish a two-level system at a temperature T 
three conditions have to be fulfilled: (i) AU ^> fcgT 
to suppress thermal hopping between the two classical 
ground states; (ii) feoi 3> ksT to observe coherent os- 
cillations; (iii) fei2 S> Seqi to have approximately a two- 
level system. For large values of AU, the energy splitting 
6eqi becomes small. Therefore, we have to find parame- 
ters where AU and &01 have reasonable values. 

From Fig. 7 we find that condition (iii) is violated if 
k is tuned too close to K c {a) whereas condition (ii) is 
violated if K becomes too large. Energy splittings Ssqi 
corresponding to approximately 25 mK look promising. 
To check condition (i) we calculate the energy difference 
AU from Eq. (47) and find values corresponding to ap- 
proximately between 100 mK and 135 mK. Therefore, we 
may expect to observe coherent quantum oscillations for 
temperatures below 25 mK which is experimentally ac- 
cessible. 

Two typical examples are shown in Fig. 8; one for a = 
1.8 and k = 0.087T (a), and one for a = 1.6 and k = 
0.687T (b), indicated by two arrows in Fig. 7. The first 
example corresponds to the case a > a c (0) while the 
second example corresponds to the case a c (7r) < a < 
a c (0). The results for the two regimes look similar: For 
the parameters of Fig. 8 (a) we find AU = 1.23 x 10~ 4 
(111 mK), Se i = 0.27 x 10~ 4 (25 mK) and 5e 12 = 1.19 X 
10~ 4 (109 mK) while for the parameters of Fig. 8 (b) we 
find AU = 1.47 x 10~ 4 (135 mK), Se Q1 = 0.27 x 10" 4 
(25 mK) and 8e 12 = 1.34 x 10~ 4 (122 mK). From Fig. 7 
we conclude that for larger values of a the energy splitting 
feoi becomes more sensitive to k. The advantages and 
disadvantages of both regimes in terms of read-out are 
discussed in the next section. 




FIG. 8. The potential energy U, Eq. (51), as a function of 
qo for a — 1.8, re = 0.087T, (a) and a = 1.6, re = 0.687T (b). 
The horizontal lines represent the corresponding eigenvalues 
eo, ei and £2 of the Schrodinger equation (53). The dashed 
lines depict the quantum-mechanical ground states uo(go) in 
arbitrary units. As in Fig. 7 the right vertical axes shows the 
temperatures corresponding to the energies of the left axis. 
In both cases the energy splitting &01 corresponds to 25 mK. 



VII. READ-OUT 

To observe coherent quantum oscillations one has to 
distinguish between the (k, k — 2tt) state described by 
/U+(x) and the (k — 2tt,k) state described by fi-(x) of 
the molecule. Therefore, it is necessary to readout its 
state. One possibility is to readout the flux associated 
with each fractional vortex in a molecule, similar to the 
earlier experiments. 10 The magnetic fluxes measured in 
these states should be distinguishable. 

In order to see if this is possible we calculate the flux 



_ $0 Z" $0 
*± = 7T / Mi.xdx = — [M±(0) - M±(-oo 



)] (56) 



measured on the left half of the LJJ for the two differ- 
ent states fJ--(x) and fJ-+(x) and the results are shown in 
Fig. 9. 

For a < a c (0) the flux difference A$ = $ + — $^ van- 
ishes at the bifurcation points K c (a). For a > a c (0) the 
value of A$ remains finite for k — > 0, as in this limit 
the two states correspond to two integer fluxons weakly 
pinned at x = ±a/2. For larger values of a the flux differ- 
ence A$ increases. Therefore, in junctions with larger a 
the two states are easier to distinguish. In this case, how- 
ever, the qubit is harder to control because the relevant 
range of k becomes smaller, see Fig. 7. 




-0.4 



FIG. 9. The fluxes $+ and defined by Eq. (56) as a 
function of re for different values of a. For each value of a the 
upper branch corresponds to $+ and the lower one to 



VIII. CONCLUSIONS 

We have presented a concept of a qubit based on a 
molecule consisting of two fractional Josephson vortices. 
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Two degenerate ground states are separated by an en- 
ergy barrier, which can be tuned during the experiment 
by changing simultaneously the values of the discontinu- 
ities. The concept may work in both linear and annular 
geometries. 

In particular, we have obtained analytical solutions for 
the stationary phases and their energies in unbiased LJ Js 
with two discontinuities. Furthermore, we have analyzed 
the stability of the stationary solutions by calculating 
the corresponding eigenmodes numerically. We have used 
these eigenmodes to map the low-energy dynamics of the 
system to the dynamics of a particle in a one-dimensional 
double-well potential and have solved the corresponding 
Schrodinger equation. 

Finally we have shown how the energy barrier can be 
tuned with the help of injector currents to reach the 
quantum limit. Our results indicate that for typical pa- 
rameters a quantum-mechanical two-level system (qubit) 
can be established for temperatures below 25 mK which 
is at the limit of modern dilution refrigerators. Ex- 
periments with such fractional vortex molecules are in 
progress in the Tubingen group. 
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Appendix: Elliptic integrals and elliptic functions 

In this appendix we introduce the elliptic integral of 
first kind and the Jacobi amplitude function and provide 
a collection of formulas which are used in the present 
paper. These relations are in accordance with Ref. 21. 



1. Elliptic integrals 

The elliptic integral of first kind is defined by 

d6 



VI - fc 2 sin 2 6 



(A.l) 



where the argument fc is called modulus. 

For <j> = 7r/ 2 it reduces to the complete elliptic integral 
of first kind 

K(fc) = F(7r/2,fc). (A.2) 
The elliptic integral F(</>, fc) obeys the symmetry relation 
F(nir±(j),k) = 2nK{k)±F(^,k), (A.3) 
where n is an integer. 



2. Jacobi amplitude function 

For k < 1 the elliptic integral (A.l) is a monotonously 
increasing function of <fi. The Jacobi amplitude function 
is the inverse of the elliptic integral. For a given value 
u = F((j), k) we have to find the corresponding integration 
limit <p. Then the Jacobi amplitude function is defined 
by 

am(u, fc) = </>. (A. 4) 

The relation 

am(u, fc) = arcsinjfc -1 sin[am(fcu, fc -1 )]} (A. 5) 

extends this definition for fc > 1. 

The Jacobi amplitude function is monotonously in- 
creasing for fc < 1 and periodic for fc > 1. Furthermore, 
it solves the differential equation 



du 



am(a, fc) 



1 — fc 2 sin 2 [ am(u, fc)] 



and for k — > 1 it reduces to 

am(«, 1) = 2arctan(e u ) — n/2. 
The three Jacobi elliptic functions 

sn(u, fc) = sin[am(w, fc)], 

and 

cn(u, fc) = cos[am(it, fc)], 

as well as 



dn(w, fc) = — [am(u, fc)], 
du 

together with the Jacobi cpsilon function 



£(u,k) = j dt dn (i,fc) 
o 



(A.6) 

(A.7) 
(A.8) 
(A.9) 

(A.10) 
(All) 



are defined in terms of the Jacobi amplitude function 
am(u, fc). 

These functions satisfy the relations 



and 



sn(w, fc) = fc -1 sn(fcit,fc -1 ), (A.12) 
cn(u,fc) = dn(fcu,fc _1 ), (A.13) 



as well as 

dn(w,fc) = cn(fcu,fc _1 ), (A.14) 

and the addition theorem 

£(ui + U2,k) = £(ui,k) + £{u 2 , fc) 

— fc 2 sn(iii, fc) sn(it2, fc) sn(ui + U2, k). (A. 15) 
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